# -*- coding: utf-8 -*-
import numpy as np
from algorithm.source_terms.source_base_formula import *
from algorithm.source_terms.flow_patterns.source_base import BaseSource

class SlugFlowCalculator(BaseSource):
    def __init__(self, standard_input):
        self.standard_input = standard_input
    
    # def __init__(self, V_L, V_g, V_sL, V_sg, V_r, V_s, alpha_L, alpha_g, rho_L, rho_g, mu_L, mu_g, D, F_E, area, p, theta, eps, sigma):
    #     self.V_L = V_L
    #     self.V_g = V_g
    #     self.V_sL = V_sL
    #     self.V_sg = V_sg
    #     self.V_r = V_r
    #     self.V_s = V_s
    #     self.alpha_L = alpha_L
    #     self.alpha_g = alpha_g
    #     self.rho_L = rho_L
    #     self.rho_g = rho_g
    #     self.mu_L = mu_L
    #     self.mu_g = mu_g
    #     self.D = D
    #     self.F_E = F_E
    #     self.area = area
    #     self.p = p
    #     self.theta = theta
    #     self.eps = eps
    #     self.sigma = sigma
    
    @staticmethod
    def cal_wet_wall_fraction(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        area: float,
        D: float,
        theta: np.ndarray,
        H_lf: np.ndarray,
        sigma: np.ndarray,
    )->np.ndarray:
        """ 计算段塞流中湿壁分数

        Args:
            V_sL (np.ndarray): 液相折算速度
            V_sg (np.ndarray): 气相折算速度
            rho_L (np.ndarray): 液体密度
            rho_g (np.ndarray): 气体密度
            area (float): 管道面面积
            D (float): 管的内径
            theta (np.ndarray): 
            H_lf (np.ndarray): 液膜持液率
            sigma (np.ndarray): 表面张力

        Returns:
            np.ndarray: 计算段塞流中湿壁分数
        """ 
        if np.cos(theta).all() < 1e-10:
            #log.debug("管道倾斜角度为近似90度,不满足计算条件。")
            print("管道倾斜角度为近似90度,不满足计算条件。")
        sigma_water = 0.07
        A_F = SlugFlowCalculator.cal_A_F(H_lf, area)
        wet_wall_fraction_init = 0.624 * (A_F / area) ** 0.374
        team1 = wet_wall_fraction_init * (sigma_water / sigma) ** 0.15
        team2 = rho_g / ((rho_L - rho_g) * np.cos(theta)) * (rho_L * V_sL ** 2 * D / sigma) ** 0.25
        team3 = (V_sg ** 2 / ((1 - H_lf) ** 2 * Physics.g.value * D)) ** 0.8
        wet_wall_fraction = team1 + team2 * team3
        return wet_wall_fraction

    @staticmethod
    def cal_A_C(H_lf: np.ndarray, area: float)->np.ndarray:
        """ 计算气芯横截面积 m^2

        Args:
            H_lf (np.ndarray): 液膜持液率
            area (float): 管道面面积

        Returns:
            np.ndarray: 气芯横截面积
        """        
        A_C = (1 - H_lf) * area
        return A_C

    @staticmethod
    def cal_A_F(H_lf: np.ndarray, area: float)->np.ndarray:
        """ 计算液膜横截面积
    
        Args:
            H_lf (np.ndarray): 液膜持液率
            area (float): 管道面面积

        Returns:
            np.ndarray: 液膜横截面积
        """        
        A_F = H_lf * area
        return A_F

    @staticmethod
    def cal_S_C(D: float, wet_wall_fraction: np.ndarray)->np.ndarray:
        """ 计算气相-壁面周长

        Args:
            D (float): 管的内径
            wet_wall_fraction (np.ndarray): 段塞流中湿壁分数

        Returns:
            np.ndarray: 气相-壁面周长
        """        
        return np.pi * D * (1 - wet_wall_fraction)

    @staticmethod
    def cal_S_F(D: float, wet_wall_fraction: np.ndarray)->np.ndarray:
        """ 计算湿壁周长

        Args:
            D (float): 管的内径
            wet_wall_fraction (np.ndarray): 段塞流中湿壁分数

        Returns:
            np.ndarray: 湿壁周长
        """        
        return np.pi * D * wet_wall_fraction

    @staticmethod
    def cal_S_CD( D: float, wet_wall_fraction: np.ndarray)->np.ndarray:
        """ 湿壁部分对应的弦长

        Args:
            D (float): 管的内径
            wet_wall_fraction (np.ndarray): 段塞流中湿壁分数

        Returns:
            np.ndarray: 湿壁部分对应的弦长
        """        
        return  D * np.sin(np.pi * wet_wall_fraction)

    @staticmethod
    def cal_A_CD(D: float, wet_wall_fraction: np.ndarray)->np.ndarray:
        """ A_CD为湿壁及其弦所包含的横截面积

        Args:
            D (float): 管的内径
            wet_wall_fraction (np.ndarray): 段塞流中湿壁分数

        Returns:
            np.ndarray: A_CD为湿壁及其弦所包含的横截面积
        """        
        A_CD = D ** 2 / 4 * (np.pi * wet_wall_fraction - np.sin(2 * np.pi * wet_wall_fraction) / 2)
        return A_CD

    @staticmethod
    def cal_S_i(D: float, wet_wall_fraction: np.ndarray, H_lf: np.ndarray, area: float)->np.ndarray:
        """ 气液界面周长

        Args:
            D (float): 管的内径
            wet_wall_fraction (np.ndarray): 段塞流中湿壁分数
            H_lf (np.ndarray): 液膜持液率
            area (float): 管道面面积

        Returns:
            np.ndarray: 气液界面周长
        """        
        S_F = SlugFlowCalculator.cal_S_F(D, wet_wall_fraction)
        A_CD = SlugFlowCalculator.cal_A_CD(D, wet_wall_fraction)
        A_F = SlugFlowCalculator.cal_A_F(H_lf, area)
        S_CD = SlugFlowCalculator.cal_S_CD(D, wet_wall_fraction)
        S_i = (S_F * (A_CD - A_F) + S_CD * A_F) / A_CD
        return S_i

    @staticmethod
    def cal_re_c(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_g: np.ndarray,
        area: float,
        D: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        H_lf: np.ndarray,
        sigma: np.ndarray,
    )->np.ndarray:
        """ 气芯的雷诺数 re_c

        Args:
            V_sL (np.ndarray): 液相折算速度
            V_sg (np.ndarray):  气相折算速度
            rho_L (np.ndarray): 液体密度
            rho_g (np.ndarray): 气体密度
            mu_g (np.ndarray): 气体黏度 气体密度
            area (float): 管道面面积
            D (float): 管的内径
            theta_centerangle (np.ndarray): 湿周圆心角
            F_E : np.ndarray 夹带率
            H_lf: np.ndarray 液膜区含液率
            sigma (np.ndarray): 表面张力

        Returns:
            np.ndarray: 气芯的雷诺数
        """        
        wet_wall_fraction = SlugFlowCalculator.cal_wet_wall_fraction(V_sL, V_sg, rho_L, rho_g, area, D, theta, H_lf, sigma)
        A_C = SlugFlowCalculator.cal_A_C(H_lf, area)
        V_c = SlugFlowCalculator.cal_V_c(V_sL, V_s, F_E, H_lf)
        S_C = SlugFlowCalculator.cal_S_C(D, wet_wall_fraction)
        S_i = SlugFlowCalculator.cal_S_i(D, wet_wall_fraction, H_lf, area)
        re_c = 4 * A_C * V_c * rho_g / ((S_C + S_i) * mu_g + 1e-15)
        return re_c

    @staticmethod
    def cal_re_F(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        area: float,
        D: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        H_lf: np.ndarray,
        sigma: np.ndarray,
    )->np.ndarray:
        """ 液膜的雷诺数 re_F

        Args:
            V_sL (np.ndarray): 液相折算速度
            V_sg (np.ndarray): 气相折算速度
            rho_L (np.ndarray): 液体密度
            rho_g (np.ndarray): 气体密度
            mu_L (np.ndarray): 液体黏度 气体密度
            area (float): 管道面面积
            D : float 管的内径 
            theta_centerangle (np.ndarray): 湿周圆心角
            F_E : np.ndarray 夹带率
            H_lf: np.ndarray 液膜区含液率
            sigma (np.ndarray): 表面张力

        Returns:
            np.ndarray: 液膜的雷诺数
        """        
        wet_wall_fraction = SlugFlowCalculator.cal_wet_wall_fraction(V_sL, V_sg, rho_L, rho_g, area, D, theta, H_lf, sigma)
        A_F = SlugFlowCalculator.cal_A_F(H_lf, area)
        V_F = SlugFlowCalculator.cal_V_F(V_sL, F_E, H_lf)
        S_F = SlugFlowCalculator.cal_S_F(D, wet_wall_fraction)
        re_F = 4 * A_F * V_F * rho_L / (S_F * mu_L + 1e-5)
        return re_F
        
    @staticmethod
    def cal_h_L(theta_centerangle: np.ndarray, D: float)-> np.ndarray:
        """ 液位高度 h_L

        Args:
            theta_centerangle (np.ndarray): 湿周圆心角
            D : float 管的内径 

        Returns:
            np.ndarray: 液位高度
        """    
        h_L = (1 - np.cos(theta_centerangle / 2)) * D /2
        return h_L

    @staticmethod
    def cal_theta_centerangle(alpha_L: np.ndarray) -> np.ndarray:
        """ " 根据液位 h_L 计算中心角 theta，公式 (1.28)

        Args:
            D : np.ndarray 管的内径
            h_L : np.ndarray 管内液位高度

        Returns:
            np.ndarray 中心角
        """
        theta_centerangle = np.pi * alpha_L + (3 * np.pi / 2) ** (1 / 3) * (1 -  2 * alpha_L + alpha_L ** (1 / 3) - (1 - alpha_L) ** (1 / 3) )
        return theta_centerangle

        # 计算液塞长度
    @staticmethod
    def cal_L_s(theta: np.ndarray, D: float) -> float:
        """ 计算液塞长度 L_s,公式为 exp(-26.8 + 28.5 * (log(D) + 1.67) ** 0.1)。

        Args:
            D : float 管的内径

        Returns:
            float 液塞长度
        """
    
        L_s = (32 * np.cos(theta) ** 2 + 16 * np.sin(theta) ** 2) * D
        return L_s

    @staticmethod
    def cal_V_t(  
        V_s: np.ndarray,
        theta: np.ndarray, 
        D: float,
    ) -> np.ndarray:
        """ 计算总速度 V_t

        Args:
            V_s : np.ndarray 混合速度 
            theta : np.ndarray 管道与水平方向的夹角
            D: float 管的内径

        Returns:
            np.ndarray 总速度
        """

        V_d = np.sqrt(Physics.g.value * D) * (0.35 * np.sin(theta) + 0.54 * np.cos(theta))
        V_t = 1.2 * V_s + V_d
        return V_t

    @staticmethod
    def cal_H_ls_estimates(
        V_sL: np.ndarray,   
        V_sg: np.ndarray,
        V_s:np.ndarray,
    ) -> np.ndarray:
        """ 计算估计值的液塞截面持液率 H_ls。

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度

        Returns:
            np.ndarray 液塞截面持液率
        """
    
        H_ls_estimates = 1 / (1 + (np.abs(V_s) / 8.66) ** 1.39)
        return H_ls_estimates

    @staticmethod
    def cal_H_lf( 
        V_sL: np.ndarray, 
        V_sg: np.ndarray,
        V_s: np.ndarray, 
        V_t: np.ndarray,
        F_E: np.ndarray,
        H_ls: np.ndarray, 
    ) -> np.ndarray:
        """ 液膜持液率 H_lf

        Args:
            H_ls (np.ndarray): 液塞持液率
            V_t (np.ndarray): 总速度
            V_s (np.ndarray): 混合速度
            V_sL (np.ndarray): 液相折算速度
            V_sg (np.ndarray): 气相折算速度
            F_E : np.ndarray 夹带率

        Returns:
            np.ndarray: 液膜持液率
        """        
        H_lf = ((H_ls * (V_t - V_s) + V_sL) * (V_sg + V_sL * F_E) - V_t * V_sL * F_E) / (V_t * V_sg + 1e-15)
        valid_indices = np.where((H_lf >= 0) & (H_lf <= 1))[0]
        if len(valid_indices) > 0:
            valid_H_lf_mean = np.mean(H_lf[valid_indices])
        else:
            valid_H_lf_mean = 0.001
            # log.debug("计算h_lf有误！")
            print("计算h_lf有误！")

        index_less_than_zero = np.where(H_lf < 0)[0]
        if len(index_less_than_zero) > 0 :
            if len(index_less_than_zero)==len(H_lf):
                H_lf[index_less_than_zero] = 0.001
            else:
                H_lf[index_less_than_zero] = valid_H_lf_mean

        index_greater_than_one = np.where(H_lf > 1)[0]
        if len(index_greater_than_one) > 0:
            H_lf[index_greater_than_one] = valid_H_lf_mean
        return H_lf

    @staticmethod
    def cal_H_lc(
        V_sL: np.ndarray,
        V_s: np.ndarray,
        F_E: np.ndarray, 
        H_lf: np.ndarray,
    )-> np.ndarray:
        """ 气芯持液率 H_lc

        Args:
            V_sL (np.ndarray): 液相折算速度
            F_E : np.ndarray 夹带率
            H_lf: np.ndarray 液膜区含液率
            V_s (np.ndarray): 混合速度

        Returns:
            np.ndarray: 气芯持液率
        """        
        H_lc = np.abs(V_sL * F_E * (1 - H_lf)) / (V_s - V_sL * (1 - F_E) + 1e-15)
        index = np.where(H_lc > 1)[0]
        if len(index) > 0:
            H_lc[index] = 1
        return H_lc

    @staticmethod
    def cal_rho_c( 
        V_s: np.ndarray, 
        V_sL: np.ndarray, 
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        F_E: np.ndarray,
        H_lf: np.ndarray, 
    ) -> np.ndarray:
        """   混合密度rho_c

        Args:
            V_s (np.ndarray): 混合速度
            V_sL (np.ndarray): 液相折算速度
            rho_L (np.ndarray): 液相密度
            rho_g (np.ndarray): 气相密度  
            F_E : np.ndarray 夹带率
            H_lf: np.ndarray 液膜区含液率
            
        Returns:
            np.ndarray: 混合密度
        """        
        H_lc = SlugFlowCalculator.cal_H_lc(V_sL, V_s, F_E, H_lf)
        rho_c = (rho_g * (1 - H_lf - H_lc) + rho_L * H_lc) / (1 - H_lf + 1e-15)
        return rho_c

    @staticmethod
    def cal_V_F(V_sL: np.ndarray, F_E: np.ndarray, H_lf: np.ndarray)-> np.ndarray:
        """ 计算液膜速度 V_F

        Args:
            V_sL np.ndarray 液相折算速度
            F_E : np.ndarray 夹带率
            H_lf: np.ndarray 液膜区含液率

        Returns:
            np.ndarray: 液膜速度
        """        
        #液膜速度
        V_F = V_sL * (1 - F_E) / (H_lf + 1e-15)
        return V_F

    @staticmethod
    def cal_lam_t(re: np.ndarray, D: float, eps: float) -> np.ndarray:
        """ 使用经验公式计算壁面摩擦因子，气液相分别和管壁的摩擦系数

        Args:
            re : np.ndarray 雷诺数
            D : float 管的内径
            eps : float 介质表面粗糙度

        Returns:
            np.ndarray 经验公式计算壁面摩擦因子
        """

        lam_t = np.zeros_like(re)

        # 雷诺数为0的情况
        mask_zero = re == 0
        lam_t[mask_zero] = 0

        # 雷诺数大于3000的情况
        mask_turbulent = re > 3000
        if np.any(mask_turbulent):
            lam_t[mask_turbulent] = SlugFlowCalculator.calculate_lam_t(re[mask_turbulent], D, eps)

        # 雷诺数小于2000的情况
        mask_laminar = re < 2000
        lam_t[mask_laminar] = 16 / re[mask_laminar]

        # 过渡区处理
        mask_transition = (re >= 2000) & (re <= 3000)
        if np.any(mask_transition):
            # 进行插值
            lam_laminar = 16 / re[mask_transition]
            lam_turbulent = SlugFlowCalculator.calculate_lam_t(re[mask_transition], D, eps)
            ratio = (re[mask_transition] - 2000) / 1000
            lam_t[mask_transition] = lam_laminar * (1 - ratio) + lam_turbulent * ratio
        return lam_t

    @staticmethod
    def calculate_lam_t(
        Re: np.ndarray, D: float, eps: float
    ) -> np.ndarray:
        """ 计算摩擦因子 lam_t, 公式 (1.23) 液相/气相与管壁

        Args:
            Re : np.ndarray 雷诺数
            d_h : np.ndarray 水力直径
            eps : np.ndarray

        Returns:
            np.ndarray 摩擦因子
        """
        # y0 = 1
        # error = 1
        # while error > 1e-6:
        #     y = 1 / (3.48 - 4 * np.log(2 * eps / D + 9.35 / Re / y0))
        #     error = np.sum(np.abs(y - y0))
        #     y0 = np.copy(y)
        # lam_t = y**2 #* 4
        lam_t = 0.0055 * (1 + ((2e4 * eps) / D + (1e6 / Re)) ** (1 / 3))
        return lam_t

    @staticmethod
    def cal_H_ls_exact_init(
        V_sL: np.ndarray, 
        V_sg: np.ndarray,
        V_s: np.ndarray,
        V_t: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        D: float,
        theta: np.ndarray,
        F_E: np.ndarray, 
        H_ls: np.ndarray, 
        sigma: np.ndarray,
        eps: float,
    )->np.ndarray:
        """ 用于计算液塞持液率HLS

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            V_t: np.ndarray  总速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            D : float 管的内径
            area : float 管道面面积
            theta : np.ndarray 管道与水平方向的夹角
            F_E: np.ndarray 夹带率
            H_ls : np.ndarray 液塞持液率
            sigma : np.ndarray 表面张力
            eps : float 介质表面粗糙度

        Returns:
            np.ndarray 液塞截面持液率
        """
        C_e = (2.5 - np.abs(np.sin(theta))) / 2
        H_lf_estimates = SlugFlowCalculator.cal_H_lf(V_sL, V_sg, V_s, V_t, F_E, H_ls)  #预估的值
        
        V_c = SlugFlowCalculator.cal_V_c(V_sL, V_s, F_E, H_lf_estimates)
        V_F = SlugFlowCalculator.cal_V_F(V_sL, F_E, H_lf_estimates)
        rho_s = rho_L * H_ls + rho_g * (1 - H_ls)
        re_s = cal_re(rho_s, D, V_s, mu_L)
        f_s = SlugFlowCalculator.cal_lam_t(re_s, D, eps)
        team1 = f_s * rho_s * V_s ** 2 / 2

        L_s = SlugFlowCalculator.cal_L_s(theta, D)
        team2 = D * rho_L * H_lf_estimates * (V_t - V_F) * (V_s - V_F) / (4 * L_s + 1e-15)

        rho_c = SlugFlowCalculator.cal_rho_c( V_s, V_sL, rho_L, rho_g, F_E, H_lf_estimates)
        V_c = SlugFlowCalculator.cal_V_c(V_sL, V_s, F_E, H_lf_estimates)
        team3 = D * rho_c * (1 - H_lf_estimates) * (V_t - V_c) * (V_s - V_c) / (4 * L_s + 1e-15)

        T_SM = (team1 + team2 + team3) / (C_e + 1e-15)

        H_ls = 1 / (1 + T_SM / (3.16 * np.sqrt((rho_L - rho_g) * Physics.g.value * sigma)))
        return H_ls

    @staticmethod
    def cal_H_ls_exact_value(
        V_sL: np.ndarray, 
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        D: float,
        theta: np.ndarray,
        F_E: np.ndarray, 
        sigma: np.ndarray,
        eps: float,
    ) -> np.ndarray:
        """ 准确液塞持液率HLS计算及液膜持液率更新

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            D : float 管的内径
            area : float 管道面面积
            theta : np.ndarray 管道与水平方向的夹角
            F_E: np.ndarray 夹带率
            sigma : np.ndarray 表面张力
            eps : float 介质表面粗糙度

        Returns:
            np.ndarray 液塞截面持液率
        """
    
        H_ls_estimates = SlugFlowCalculator.cal_H_ls_estimates(V_sL, V_sg, V_s)
        V_t = SlugFlowCalculator.cal_V_t(V_s, theta, D) 
    
        H_ls = SlugFlowCalculator.cal_H_ls_exact_init(V_sL,V_sg, V_s, V_t, rho_L, rho_g, mu_L, D, theta, F_E, H_ls_estimates, sigma, eps)
        
        index = np.where(H_ls < 0.24)[0]
        if len(index) > 0:
            H_ls[index] = 0.24
        index = np.where(H_ls > 1)[0]
        if len(index) > 0:
            if index[0]==1:
                H_ls[index] = 0.5
            else:
                H_ls[index] = H_ls[index[0]-1]
        
        H_ls = SlugFlowCalculator.cal_H_ls_exact_init(V_sL,V_sg, V_s, V_t, rho_L, rho_g, mu_L, D, theta, F_E, H_ls, sigma, eps)
        
        index = np.where(H_ls < 0.24)[0]
        if len(index) > 0:
            H_ls[index] = 0.24
        index = np.where(H_ls > 1)[0]
        if len(index) > 0:
            if index[0]==1:
                H_ls[index] = 0.5
            else:
                H_ls[index] = H_ls[index[0]-1]
        H_lf = SlugFlowCalculator.cal_H_lf(V_sL, V_sg, V_s, V_t, F_E, H_ls)
        # index = np.where(H_lf < 0)[0]
        # if len(index) > 0:
        #     H_lf[index] = 0.0001
        # index = np.where(H_lf > 1)[0]
        # if len(index) > 0:
        #     H_lf[index] = 0.78

        return H_ls, H_lf

    @staticmethod
    def cal_V_gs(
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        theta: np.ndarray,
        H_ls: np.ndarray,
        sigma: np.ndarray,  
    ) -> np.ndarray:
        """ 计算气泡区域流速 V_gs。

        Args:
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            theta : np.ndarry 管道与水平方向的夹角
            H_ls : np.ndarray 液塞持液率
            sigma : np.ndarray 表面张力

        Returns:
            np.ndarray 气泡区域流速
        """

        V_gs = 1.2 * V_s + 1.53 * (
            sigma * Physics.g.value * (rho_L - rho_g) / (rho_L ** 2 + 1e-15)
        ) ** 0.25 * H_ls ** 0.1 * np.sin(theta)
        return V_gs

    @staticmethod
    def cal_V_ls(
        V_s: np.ndarray,
        theta: np.ndarray,
        rho_L: np.ndarray,
        sigma: np.ndarray,
        rho_g: np.ndarray,
        H_ls: np.ndarray,
    ) -> np.ndarray:
        """ 该函数根据输入参数计算液塞流速度 V_ls。

        Args:
            V_s : np.ndarray 混合速度
            theta : np.ndarray 管道与水平方向的夹角
            rho_L : np.ndarray 液体密度
            sigma : np.ndarray 表面张力
            rho_g : np.ndarray 气体密度
            H_ls : np.ndarray 液塞持液率

        Returns:
            np.ndarray 液塞流速度
        """
        V_gs = SlugFlowCalculator.cal_V_gs(V_s, rho_L, rho_g, theta, H_ls, sigma)
        V_ls = V_gs + (V_s - V_gs) / (H_ls + 1e-15)
        return V_ls

    @staticmethod
    def cal_L_u(
        V_s: np.ndarray,
        V_sL: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        D: float,
        theta: np.ndarray,
        H_lf:np.ndarray, 
        H_ls: np.ndarray,
        sigma: np.ndarray,
    ) -> np.ndarray:
        """ 计算单位长度液塞长度 L_u

        Args:
            V_sL : np.ndarray 液相折算速度
            V_s: np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            H_lf : np.ndarray 液膜持液率
            H_ls : np.ndarray 液塞持液率
            sigma : np.ndarray 表面张力

        Returns:
            np.ndarray 液塞长度
        """
        
        V_ls = SlugFlowCalculator.cal_V_ls(V_s, theta, rho_L, sigma, rho_g, H_ls)
        L_s = SlugFlowCalculator.cal_L_s(theta, D)

        V_t = SlugFlowCalculator.cal_V_t(V_s, theta, D) 
        
        V_lf = V_t - (V_t - V_ls) * H_ls / (H_lf + 1e-15)
        L_u = L_s * (V_ls * H_ls - V_lf * H_lf) / (V_sL - V_lf * H_lf + 1e-15)
        return L_u

    @staticmethod
    def cal_tau_fw(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        D: float,
        area: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        H_lf: np.ndarray,
        sigma: np.ndarray,
        eps: float,
    ) -> np.ndarray:
        """ 计算液塞区域剪切应力 tau_ls

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            area : float 管道面面积
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            sigma : np.ndarray 表面张力
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            eps : float 介质表面粗糙度

        Returns:
            np.ndarray 液塞区域剪切应力
        """
        f_fw = SlugFlowCalculator.cal_slug_flow_wall_fric(
            V_sL,
            V_sg, 
            V_s,
            rho_L,
            rho_g,
            mu_L,
            mu_g,
            area,
            D,
            theta,
            F_E,
            H_lf,
            sigma,
            eps,
            flag="liquid",
        )
        
        U_F = V_sL * (1 - F_E) / (H_lf + 1e-15)
        tau_fw = 0.5 * f_fw * rho_L * U_F ** 2
        return tau_fw

    @staticmethod
    def cal_tau_cw(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        D: float,
        area: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        H_lf: np.ndarray,
        sigma: np.ndarray, 
        eps: float,
    ) -> np.ndarray:
        """ 计算气泡区域剪切应力 tau_cw

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            area : float 管道面面积
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            sigma : np.ndarray 表面张力
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            eps : float 介质表面粗糙度

        Returns:
            np.float64: tau_cw numpy
        """
        f_cw = SlugFlowCalculator.cal_slug_flow_wall_fric(
            V_sL,
            V_sg, 
            V_s,
            rho_L,
            rho_g,
            mu_L,
            mu_g,
            area,
            D,
            theta,
            F_E,
            H_lf,
            sigma,
            eps,
            flag = "gas",
        )
        
        V_c = SlugFlowCalculator.cal_V_c(V_sL, V_s, F_E, H_lf)
        tau_cw = 0.5 * f_cw * rho_g * V_c ** 2
        # log.debug("函数V_gs：%s",V_gs)
        return tau_cw

    @staticmethod
    def cal_V_c(
        V_sL: np.ndarray, 
        V_s: np.ndarray,  
        F_E: np.ndarray,    
        H_lf: np.ndarray
    )-> np.ndarray:
        """ 计算气芯速度V_c

        Args:
            V_sL (np.ndarray): 液相折算速度       
            F_E : np.ndarray 夹带率
            H_lf: np.ndarray 液膜区含液率

        Returns:
            np.ndarray: 气芯速度V_c
        """        
        V_c = (V_s - V_sL * (1 - F_E)) / (1 - H_lf + 1e-15)
        return V_c

    @staticmethod
    def cal_slug_flow_wall_fric(
        V_sL: np.ndarray,
        V_sg: np.ndarray, 
        V_s: np.ndarray, 
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        area: float,
        D: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        H_lf: np.ndarray,
        sigma: np.ndarray,
        eps: float,
        flag="liquid",
    )-> np.ndarray:
        """ 计算摩擦力系数

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            area : float 管道面面积
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            sigma : np.ndarray 表面张力
            eps : float 介质表面粗糙度
            flag (str, optional): _description_. Defaults to "liquid".

        Returns:
            计算摩擦力系数 lam_t
        """    
        if flag == "liquid":
            re = SlugFlowCalculator.cal_re_F(
                V_sL,
                V_sg, 
                rho_L,
                rho_g,
                mu_L,
                area,
                D,
                theta,
                F_E,
                H_lf,
                sigma)
        elif flag == "gas":
            re = SlugFlowCalculator.cal_re_c(
                    V_sL,
                    V_sg, 
                    V_s,  
                    rho_L,
                    rho_g,
                    mu_g,
                    area,
                    D,
                    theta,
                    F_E,
                    H_lf,
                    sigma)
        lam_t = SlugFlowCalculator.cal_lam_t(re, D, eps)    
        return lam_t

    @staticmethod
    def cal_slug_flow_border_fric( 
        D: np.ndarray,
        V_sg: np.ndarray,
        p: np.ndarray,
        h_lf: np.ndarray,
        lam_g: np.ndarray,
    ) -> np.ndarray:
        """ 求气液相之间的界面摩擦系数

        Args:
        
            D : float 管的内径  
            V_sg : np.ndarray 气相折算速度
            p : np.ndarray 压强
            h_L : np.ndarray 管内液位高度
            lam_g : np.ndarray  气壁摩阻系数

        Returns:
            np.ndarray 界面摩擦系数
        """
        # 保持原有注释
        v_sgt = 5 * (101325 / p) ** 0.5  # 使用矢量化运算计算v_sgt

        #V_sg = V_g * alpha_g  # 计算V_sg
        condition = V_sg <= v_sgt  # 计算条件

        lam_i = np.where(
            condition,
            lam_g,
            lam_g * (1 + 14.3 * (h_lf ** 0.5) * (V_sg / v_sgt - 1)),
        )  # 使用np.where替代for循环
        # lam_i = 0.02*(1+75*(1 - V_g)) # TODO:新公式
        return lam_i

    @staticmethod
    def cal_tau_i(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        D: float,
        area: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        p: np.ndarray,
        h_L: np.ndarray,
        H_lf: np.ndarray,
        sigma: np.ndarray,
        eps: float,
    ) -> np.ndarray:
        """ 计算液膜剪切应力 tau_gs

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            area : float 管道面面积
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            p: np.ndarray 压强
            h_L: np.ndarray 液膜区液位高度
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            sigma : np.ndarray 表面张力
            eps : float 介质表面粗糙度

        Returns:
            np.float64: tau_gs numpy
        """
    
        f_cw = SlugFlowCalculator.cal_slug_flow_wall_fric(
            V_sL,
            V_sg, 
            V_s,
            rho_L,
            rho_g,
            mu_L,
            mu_g,
            area,
            D,
            theta,
            F_E,
            H_lf,
            sigma,
            eps,
            flag="liquid",
            )

        f_i = SlugFlowCalculator.cal_slug_flow_border_fric(
                D,
                V_sg,
                p,
                H_lf,#h_L, #H_lf
                f_cw,
            )
        V_F = SlugFlowCalculator.cal_V_F(V_sL, F_E, H_lf)    
        V_c = SlugFlowCalculator.cal_V_c(V_sL, V_s, F_E, H_lf)
        rho_c = SlugFlowCalculator.cal_rho_c(V_s, V_sL, rho_L, rho_g, F_E, H_lf)
        tau_i = 0.5 * f_i * rho_c * (V_c - V_F) ** 2
        #tau_i = 0.5 * f_i * rho_g * (V_c - V_F) ** 2
        # log.debug("函数V_gs：%s",V_gs)
        return tau_i

    @staticmethod
    def cal_tau_s(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s:np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        D: float,
        area: float,
        theta: np.ndarray, 
        F_E: np.ndarray,
        H_lf: np.ndarray,
        H_ls: np.ndarray,
        sigma: np.ndarray,
        eps: float,
    ) -> np.ndarray:
        """ 计算液塞区域剪切应力 tau_ls

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            area : float 管道面面积
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            sigma : np.ndarray 表面张力
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            eps : float 介质表面粗糙度

        Returns:
            np.ndarray 液塞区域剪切应力
        """
        
        f_s = SlugFlowCalculator.cal_slug_flow_wall_fric(
            V_sL,
            V_sg, 
            V_s,
            rho_L,
            rho_g,
            mu_L,
            mu_g,
            area,
            D,
            theta,
            F_E,
            H_lf,
            sigma,
            eps,
            flag="liquid",
        )
    
        rho_s = rho_L *  H_ls + rho_g * (1 - H_ls)
        tau_s = 0.5 * f_s * rho_s * V_s ** 2
        return tau_s

    @staticmethod
    def cal_slug_flow_gamma_gw(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        D: float,
        area: float,
        theta: np.ndarray, 
        F_E: np.ndarray,
        H_lf: np.ndarray,
        H_ls: np.ndarray, 
        sigma: np.ndarray,
        eps: float,
    ) -> np.ndarray:
        """  gamma_gw (np.ndarray float)

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            area : float 管道面面积
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            sigma : np.ndarray 表面张力
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            eps : float 介质表面粗糙度

        Returns:
            np.ndarray gamma_gw
        """
        
        tau_cw = SlugFlowCalculator.cal_tau_cw(
            V_sL, V_sg, V_s, rho_L, rho_g,  mu_L, mu_g, D, area, theta, F_E, H_lf, sigma, eps
        )

        wet_wall_fraction = SlugFlowCalculator.cal_wet_wall_fraction(V_sL, V_sg, rho_L, rho_g, area, D, theta, H_lf, sigma)
        S_C = SlugFlowCalculator.cal_S_C(D, wet_wall_fraction)
        L_s = SlugFlowCalculator.cal_L_s(theta, D)

        L_u = SlugFlowCalculator.cal_L_u(
            V_s, V_sL, rho_L, rho_g, D, theta, H_lf, H_ls, sigma
        )
        L_f = L_u - L_s
        gamma_gw =  -L_f / L_u * tau_cw * S_C
        
        # log.debug("函数中S_gs：%s",S_gs)
        # log.debug("函数中tau_gs：%s",tau_gs)
        # log.debug("函数中V_gf：%s",V_gf)
        # log.debug("函数中tau_gf：%s",tau_gf)
        # log.debug("函数中S_gf：%s",S_gf)
        return gamma_gw

    @staticmethod
    def cal_slug_flow_gamma_Lw(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        D: float,
        area: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        H_lf: np.ndarray,
        H_ls: np.ndarray,
        sigma: np.ndarray,
        eps: float,
    ) -> np.ndarray:
        """  计算 gamma_Lw

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            area : float 管道面面积
            D : float 管的内径
            theta : np.ndarray 管道与水平方向的夹角
            sigma : np.ndarray 表面张力
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            eps : float 介质表面粗糙度

        Returns:
            np.ndarray gamma_Lw
        """
        tau_fw = SlugFlowCalculator.cal_tau_fw(
            V_sL, V_sg, V_s, rho_L, rho_g, mu_L, mu_g, D, area, theta, F_E, H_lf, sigma, eps)

        tau_s = SlugFlowCalculator.cal_tau_s(
            V_sL, V_sg, V_s, rho_L, rho_g, mu_L, mu_g, D, area, theta, F_E, H_lf, H_ls, sigma, eps)

        L_s = SlugFlowCalculator.cal_L_s(theta, D)
        L_u = SlugFlowCalculator.cal_L_u(
            V_s, V_sL, rho_L, rho_g, D, theta, H_lf, H_ls, sigma
        )
    
        L_f = L_u - L_s

        wet_wall_fraction = SlugFlowCalculator.cal_wet_wall_fraction(V_sL, V_sg, rho_L, rho_g, area, D, theta, H_lf, sigma)
        
        S_F = SlugFlowCalculator.cal_S_F(D, wet_wall_fraction)
        gamma_Lw = L_s / L_u * tau_s * np.pi * D + L_f / L_u * tau_fw * S_F

        return gamma_Lw

    @staticmethod
    def cal_slug_flow_gamma_i(
        V_sL: np.ndarray,
        V_sg: np.ndarray,
        V_s: np.ndarray,
        rho_L: np.ndarray,
        rho_g: np.ndarray,
        mu_L: np.ndarray,
        mu_g: np.ndarray,
        D: float,
        area: float,
        theta: np.ndarray,
        F_E: np.ndarray,
        p: np.ndarray,
        h_L: np.ndarray,
        H_lf: np.ndarray,
        H_ls: np.ndarray,
        sigma: np.ndarray,
        eps: float,
    ) -> np.ndarray:
        """ 计算 gamma_i

        Args:
            V_sL : np.ndarray 液相折算速度
            V_sg : np.ndarray 气相折算速度
            V_s : np.ndarray 混合速度
            alpha_g : np.ndarray 管道中连续的气体的截面分数
            rho_L : np.ndarray 液体密度
            rho_g : np.ndarray 气体密度
            mu_L : np.ndarray  液体黏度
            mu_g : np.ndarray  气体黏度
            D : float 管的内径
            area : float 管道面面积 
            theta : np.ndarray 管道与水平方向的夹角
            p: np.ndarray 压强
            h_L: np.ndarray 液膜区液位高度
            F_E: np.ndarray 夹带率
            H_ls: np.ndarray 液塞截面持液率
            H_lf: np.ndarray 液膜区含液率
            sigma : np.ndarray 表面张力
            eps : float 介质表面粗糙度

        Returns:
            np.float64: gamma_i numpy
        """ 
    
        tau_i = SlugFlowCalculator.cal_tau_i(
            V_sL, V_sg, V_s, rho_L, rho_g, mu_L,mu_g, D, area, theta, F_E, p, h_L, H_lf, sigma, eps
        )
        
        L_s = SlugFlowCalculator.cal_L_s(theta, D)

        L_u = SlugFlowCalculator.cal_L_u(
            V_s, V_sL, rho_L, rho_g, D, theta, H_lf, H_ls, sigma
        )
        L_f = L_u - L_s
        
        #S_i = cal_slug_flow_S_i(D, L_f, L_u)
        wet_wall_fraction = SlugFlowCalculator.cal_wet_wall_fraction(V_sL, V_sg, rho_L, rho_g, area, D, theta, H_lf, sigma)
        S_i = SlugFlowCalculator.cal_S_i(D, wet_wall_fraction, H_lf, area)
        gamma_i = L_f / L_u * tau_i * S_i
        #gamma_i =  tau_i * S_i
        # log.debug("函数中C_D:%s",C_D)
        return gamma_i


    def calculate(self):
        self.theta_centerangle = SlugFlowCalculator.cal_theta_centerangle(
                self.alpha_L
            )

        self.h_L = SlugFlowCalculator.cal_h_L(
            self.theta_centerangle, self.D
        )

        self.H_ls, self.H_lf = SlugFlowCalculator.cal_H_ls_exact_value( 
            self.V_sL, 
            self.V_sg, 
            self.V_s, 
            self.rho_L, 
            self.rho_g, 
            self.mu_L, 
            self.D, 
            self.theta,  
            self.F_E, 
            self.sigma, 
            self.eps,
            )
        
        self.gamma_gw = SlugFlowCalculator.cal_slug_flow_gamma_gw(    
            self.V_sL,
            self.V_sg,
            self.V_s,
            self.rho_L,
            self.rho_g,
            self.mu_L,
            self.mu_g,
            self.D,
            self.area,
            self.theta, 
            self.F_E, 
            self.H_lf,  
            self.H_ls, 
            self.sigma,  
            self.eps,
        )

        self.gamma_Lw = SlugFlowCalculator.cal_slug_flow_gamma_Lw(
            self.V_sL,
            self.V_sg,
            self.V_s,
            self.rho_L,
            self.rho_g,
            self.mu_L,
            self.mu_g,
            self.D,
            self.area,
            self.theta,
            self.F_E,  
            self.H_lf,
            self.H_ls, 
            self.sigma, 
            self.eps,
        )

        self.gamma_i = SlugFlowCalculator.cal_slug_flow_gamma_i(
            self.V_sL,
            self.V_sg,
            self.V_s,
            self.rho_L,
            self.rho_g,
            self.mu_L,
            self.mu_g,
            self.D,
            self.area,
            self.theta,
            self.F_E,  
            self.p,
            self.h_L,
            self.H_lf,  
            self.H_ls,
            self.sigma,
            self.eps,
        )

        return self.gamma_gw, self.gamma_Lw, self.gamma_i